module complete_elliptic
!
!   Complete elliptic integral of the first kind
!   Input:  mc = 1-m  "Complementary parameter m"   0 <= mc <= 1
!   Output: elk
!   Dependencies: None
!
!   Please note that the parameter m is related to the parameter k used in other notation as
!           m = k**2
!
!   This implementation written by Magnus Lundmark as a conversion and 
!   extension of Prof. Fukushimas original F77 code.
!
!   Reference: T. Fukushima, (2009) Celest. Mech. Dyn. Astron. 105, 305-328
!   "Fast Computation of Complete Elliptic Integrals and Jacobian Elliptic Functions"
!   Author: T. Fukushima Toshio.Fukushima@nao.ac.jp
!
!   Copyright 2012 Magnus Lundmark
  use select_precision
  use constants
  implicit none
  
  interface CompleteEllipticFirstKind
    module procedure sp_elk , dp_elk , tp_elk , qp_elk
  end interface CompleteEllipticFirstKind
  
contains  
  pure function sp_elk( mc )
    !Functions variables
    real(sp),intent(in)   :: mc
    real(sp)                  :: sp_elk
    !Working variables
    real(sp)                  :: m , mx , P , Q , kkc , nome 
    !Parameters
    real(sp),parameter :: D01 =     real(1_i8,sp)/real(16_i8,sp)
    real(sp),parameter :: D02 =     real(1_i8,sp)/real(32_i8,sp) 
    real(sp),parameter :: D03 =    real(21_i8,sp)/real(1024_i8,sp)
    real(sp),parameter :: D04 =    real(31_i8,sp)/real(2048_i8,sp) 
    real(sp),parameter :: D05 =  real(6257_i8,sp)/real(524288_i8,sp)
    real(sp),parameter :: D06 = real(10293_i8,sp)/real(1048576_i8,sp) 

          m = 1.0_sp - mc
argument: if( mc < 1.05e-8_sp) then
              sp_elk = LN4 - 0.5_sp*log(mc)
          elseif(mc < 0.1_sp) then
              nome = mc*(D01 + mc*(D02 + mc*(D03 + mc*(D04 + mc*(D05 + mc*D06)))))
              mx = mc-0.05_sp
             !   K'
              kkc=1.59100345_sp + mx*(&
                  0.41600074_sp + mx*(&
                  0.24579151_sp + mx*(&
                  0.17948148_sp + mx*(&
                  0.14455606_sp))))
              sp_elk = -kkc*PIINV*log(nome)
    elseif(m < 0.1_sp) then
              mx = m - 0.05_sp
              sp_elk=1.59100345_sp + mx*(&
                     0.41600074_sp + mx*(&
                     0.24579151_sp + mx*(&
                     0.17948148_sp + mx*(&
                     0.14455606_sp))))
    elseif(m < 0.2_sp) then
              mx=m-0.15_sp
              sp_elk=1.63525673_sp + mx*(&
                     0.47119063_sp + mx*(&
                     0.30972841_sp + mx*(&
                     0.25220831_sp + mx*(&
                     0.22672562_sp))))
    elseif(m < 0.3_sp) then
              mx=m-0.25_sp
              sp_elk=1.68575035_sp + mx*(&
                     0.54173185_sp + mx*(&
                     0.40152444_sp + mx*(&
                     0.36964247_sp + mx*(&
                     0.37606072_sp))))
    elseif(m < 0.4_sp) then
              mx=m-0.35_sp
              sp_elk=1.74435060_sp + mx*(&
                     0.63486428_sp + mx*(&
                     0.53984256_sp + mx*(&
                     0.57189271_sp + mx*(&
                     0.67029514_sp + mx*(&
                     0.83258659_sp)))))
    elseif(m < 0.5_sp) then
              mx=m-0.45_sp
              sp_elk=1.81388394_sp + mx*(&
                     0.76316325_sp + mx*(&
                     0.76192861_sp + mx*(&
                     0.95107465_sp + mx*(&
                     1.31518067_sp + mx*(&
                     1.92856069_sp)))))
    elseif(m < 0.6_sp) then
              mx=m-0.55_sp
              sp_elk=1.89892491_sp + mx*(&
                     0.95052179_sp + mx*(&
                     1.15107759_sp + mx*(&
                     1.75023911_sp + mx*(&
                     2.95267681_sp + mx*(&
                     5.28580040_sp)))))
    elseif(m < 0.7_sp) then
              mx=m-0.65_sp
              sp_elk=2.00759840_sp + mx*(&
                     1.24845723_sp + mx*(&
                     1.92623466_sp + mx*(&
                     3.75128964_sp + mx*(&
                     8.11994455_sp + mx*(&
                     18.6657213_sp + mx*(&
                     44.6039248_sp))))))
    elseif(m < 0.8_sp) then
              mx=m-0.75_sp
              sp_elk=2.15651565_sp + mx*(&
                     1.79180564_sp + mx*(&
                     3.82675129_sp + mx*(&
                     10.3867247_sp + mx*(&
                     31.4033141_sp + mx*(&
                     100.923704_sp + mx*(&
                     337.326828_sp + mx*(&
                     1158.70793_sp)))))))
    elseif(m < 0.85_sp) then
              mx=m-0.825_sp
              sp_elk=2.31812262_sp + mx*(&
                     2.61692015_sp + mx*(&
                     7.89793508_sp + mx*(&
                     30.5023972_sp + mx*(&
                     131.486937_sp + mx*(&
                     602.984764_sp + mx*(&
                     2877.02462_sp))))))
      else
              mx=m-0.875_sp
              sp_elk=2.47359617_sp + mx*(&
                     3.72762424_sp + mx*(&
                     15.6073930_sp + mx*(&
                     84.1285084_sp + mx*(&
                     506.981820_sp + mx*(&
                     3252.27706_sp + mx*(&
                     21713.2424_sp + mx*(&
                     149037.045_sp)))))))
    endif argument
    return
  end function sp_elk
  
  
  pure function tp_elk( mc )
    real(tp),intent(in)      :: mc
    real(tp)                 :: tp_elk
    !Working variables
    real(tp)                 :: m , mx , P , Q , kkc , nome 
    !Parameters
    real(tp),parameter       :: very_tiny = 1.e-99_tp
    real(tp),parameter :: D01 =                      real(1_i8,tp)/real(16_i8,tp)
    real(tp),parameter :: D02 =                      real(1_i8,tp)/real(32_i8,tp) 
    real(tp),parameter :: D03 =                     real(21_i8,tp)/real(1024_i8,tp)
    real(tp),parameter :: D04 =                     real(31_i8,tp)/real(2048_i8,tp) 
    real(tp),parameter :: D05 =                   real(6257_i8,tp)/real(524288_i8,tp)
    real(tp),parameter :: D06 =                  real(10293_i8,tp)/real(1048576_i8,tp) 
    real(tp),parameter :: D07 =                 real(279025_i8,tp)/real(33554432_i8,tp)
    real(tp),parameter :: D08 =                 real(483127_i8,tp)/real(67108864_i8,tp)
    real(tp),parameter :: D09 =              real(435506703_i8,tp)/real(68719476736_i8,tp)
    real(tp),parameter :: D10 =              real(776957575_i8,tp)/real(137438953472_i8,tp)
    real(tp),parameter :: D11 =            real(22417045555_i8,tp)/real(4398046511104_i8,tp)
    real(tp),parameter :: D12 =            real(40784671953_i8,tp)/real(8796093022208_i8,tp)
    real(tp),parameter :: D13 =          real(9569130097211_i8,tp)/real(2251799813685248_i8,tp)
    real(tp),parameter :: D14 =         real(17652604545791_i8,tp)/real(4503599627370496_i8,tp)
    real(tp),parameter :: D15 =        real(523910972020563_i8,tp)/real(144115188075855872_i8,tp)
    real(tp),parameter :: D16 =        real(976501268709949_i8,tp)/real(288230376151711744_i8,tp)
    !!!Integer overflow for integer_kind=8 for the following coefficients
    real(tp),parameter :: D17 =             935823746406530603._tp/295147905179352825856._tp
    real(tp),parameter :: D18 =            1758220447807291611._tp/590295810358705651712._tp
    !Not needed for 80-bit accuracy
    !real(tp),parameter :: D19 =           53030538453624441751._tp/18889465931478580854784._tp
    !real(tp),parameter :: D20 =          100268465197007602645._tp/37778931862957161709568._tp
         
    m = 1.0_tp - mc
arg:    if ( abs(m) < 1.e-16_tp ) then
          tp_elk = PIHALF
        elseif ( mc < very_tiny ) then
          !tp_elk = LN4 - 0.5_tp*log( real(very_tiny,dp) )
          tp_elk = real(LN4,tp) - 0.5_tp*log( very_tiny )
        elseif ( mc < 1.11e-16_tp ) then
          tp_elk = LN4 - 0.5_tp*log( real(mc,dp) )
          !tp_elk = real(LN4,tp) - 0.5_tp*log( mc )
        elseif ( mc < 0.1_tp ) then
          nome = mc*(D01 + mc*(D02 + mc*(D03 + mc*(D04 + mc*(D05 + mc*(D06 &
               + mc*(D07 + mc*(D08 + mc*(D09 + mc*(D10 + mc*(D11 + mc*(D12 &
               + mc*(D13 + mc*(D14 + mc*(D15 + mc*(D16 + mc*(D17 + mc*D18)))))))))))))))))
          mx   = mc - 0.05_tp
!
!      K'
!
        kkc=1.591003453790792180_tp+mx*( &
            0.416000743991786912_tp+mx*( &
            0.245791514264103415_tp+mx*( &
            0.179481482914906162_tp+mx*( &
            0.144556057087555150_tp+mx*( &
            0.123200993312427711_tp+mx*( &
            0.108938811574293531_tp+mx*( &
            0.098853409871592910_tp+mx*( &
            0.091439629201749751_tp+mx*( &
            0.085842591595413900_tp+mx*( &
            0.081541118718303215_tp))))))))))
        tp_elk = -kkc*PIINV*log( real(nome,dp) )
        !tp_elk = -kkc*real(PIINV,tp)*log( nome )
    elseif(m <= 0.1_tp) then
            mx =  m-0.05_tp
            tp_elk = 1.591003453790792180_tp+mx*( &
                  0.416000743991786912_tp+mx*( &
                  0.245791514264103415_tp+mx*( &
                  0.179481482914906162_tp+mx*( &
                  0.144556057087555150_tp+mx*( &
                  0.123200993312427711_tp+mx*( &
                  0.108938811574293531_tp+mx*( &
                  0.098853409871592910_tp+mx*( &
                  0.091439629201749751_tp+mx*( &
                  0.085842591595413900_tp+mx*( &
                  0.081541118718303215_tp))))))))))
    elseif(m <= 0.2_tp) then
            mx  = m-0.15_tp
            tp_elk = 1.635256732264579992_tp+mx*( &
                  0.471190626148732291_tp+mx*( &
                  0.309728410831499587_tp+mx*( &
                  0.252208311773135699_tp+mx*( &
                  0.226725623219684650_tp+mx*( &
                  0.215774446729585976_tp+mx*( &
                  0.213108771877348910_tp+mx*( &
                  0.216029124605188282_tp+mx*( &
                  0.223255831633057896_tp+mx*( &
                  0.234180501294209925_tp+mx*( &
                  0.248557682972264071_tp+mx*( &
                  0.266363809892617521_tp+mx*( &
                  0.287728452156114668_tp))))))))))))
    elseif(m <= 0.3_tp) then
            mx =  m-0.25_tp
            tp_elk = 1.685750354812596043_tp+mx*( &
                  0.541731848613280329_tp+mx*( &
                  0.401524438390690257_tp+mx*( &
                  0.369642473420889090_tp+mx*( &
                  0.376060715354583645_tp+mx*( &
                  0.405235887085125919_tp+mx*( &
                  0.453294381753999079_tp+mx*( &
                  0.520518947651184205_tp+mx*( &
                  0.609426039204995055_tp+mx*( &
                  0.724263522282908870_tp+mx*( &
                  0.871013847709812357_tp+mx*( &
                  1.057652872753547036_tp)))))))))))
    elseif(m <= 0.4_tp) then
            mx =  m-0.35_tp
            tp_elk = 1.744350597225613243_tp+mx*( &
                  0.634864275371935304_tp+mx*( &
                  0.539842564164445538_tp+mx*( &
                  0.571892705193787391_tp+mx*( &
                  0.670295136265406100_tp+mx*( &
                  0.832586590010977199_tp+mx*( &
                  1.073857448247933265_tp+mx*( &
                  1.422091460675497751_tp+mx*( &
                  1.920387183402304829_tp+mx*( &
                  2.632552548331654201_tp+mx*( &
                  3.652109747319039160_tp+mx*( &
                  5.115867135558865806_tp+mx*( &
                  7.224080007363877411_tp))))))))))))
    elseif(m <= 0.5_tp) then
            mx  = m-0.45_tp
            tp_elk = 1.813883936816982644_tp+mx*( &
                  0.763163245700557246_tp+mx*( &
                  0.761928605321595831_tp+mx*( &
                  0.951074653668427927_tp+mx*( &
                  1.315180671703161215_tp+mx*( &
                  1.928560693477410941_tp+mx*( &
                  2.937509342531378755_tp+mx*( &
                  4.594894405442878062_tp+mx*( &
                  7.330071221881720772_tp+mx*( &
                  11.87151259742530180_tp+mx*( &
                  19.45851374822937738_tp+mx*( &
                  32.20638657246426863_tp+mx*( &
                  53.73749198700554656_tp+mx*( &
                  90.27388602940998849_tp)))))))))))))
    elseif(m <= 0.6_tp) then
            mx  = m-0.55_tp
            tp_elk = 1.898924910271553526_tp+mx*( &
                  0.950521794618244435_tp+mx*( &
                  1.151077589959015808_tp+mx*( &
                  1.750239106986300540_tp+mx*( &
                  2.952676812636875180_tp+mx*( &
                  5.285800396121450889_tp+mx*( &
                  9.832485716659979747_tp+mx*( &
                  18.78714868327559562_tp+mx*( &
                  36.61468615273698145_tp+mx*( &
                  72.45292395127771801_tp+mx*( &
                  145.1079577347069102_tp+mx*( &
                  293.4786396308497026_tp+mx*( &
                  598.3851815055010179_tp+mx*( &
                  1228.420013075863451_tp+mx*( &
                  2536.529755382764488_tp))))))))))))))
    elseif(m <= 0.7_tp) then
            mx =  m-0.65_tp
            tp_elk = 2.007598398424376302_tp+mx*( &
                  1.248457231212347337_tp+mx*( &
                  1.926234657076479729_tp+mx*( &
                  3.751289640087587680_tp+mx*( &
                  8.119944554932045802_tp+mx*( &
                  18.66572130873555361_tp+mx*( &
                  44.60392484291437063_tp+mx*( &
                  109.5092054309498377_tp+mx*( &
                  274.2779548232413480_tp+mx*( &
                  697.5598008606326163_tp+mx*( &
                  1795.716014500247129_tp+mx*( &
                  4668.381716790389910_tp+mx*( &
                  12235.76246813664335_tp+mx*( &
                  32290.17809718320818_tp+mx*( &
                  85713.07608195964685_tp+mx*( &
                  228672.1890493117096_tp+mx*( &
                  612757.2711915852774_tp))))))))))))))))
    elseif(m <= 0.8_tp) then
            mx =  m-0.75_tp
            tp_elk = 2.156515647499643235_tp+mx*( &
                  1.791805641849463243_tp+mx*( &
                  3.826751287465713147_tp+mx*( &
                  10.38672468363797208_tp+mx*( &
                  31.40331405468070290_tp+mx*( &
                  100.9237039498695416_tp+mx*( &
                  337.3268282632272897_tp+mx*( &
                  1158.707930567827917_tp+mx*( &
                  4060.990742193632092_tp+mx*( &
                  14454.00184034344795_tp+mx*( &
                  52076.66107599404803_tp+mx*( &
                  189493.6591462156887_tp+mx*( &
                  695184.5762413896145_tp+mx*( &
                  2.567994048255284686d6+mx*( &
                  9.541921966748386322d6+mx*( &
                  3.563492744218076174d7+mx*( &
                  1.336692984612040871d8+mx*( &
                  5.033521866866284541d8+mx*( &
                  1.901975729538660119d9+mx*( &
                  7.208915015330103756d9)))))))))))))))))))
    elseif(m <= 0.85_tp) then
            mx =  m-0.825_tp
            tp_elk=  2.318122621712510589_tp+mx*( &
                  2.616920150291232841_tp+mx*( &
                  7.897935075731355823_tp+mx*( &
                  30.50239715446672327_tp+mx*( &
                  131.4869365523528456_tp+mx*( &
                  602.9847637356491617_tp+mx*( &
                  2877.024617809972641_tp+mx*( &
                  14110.51991915180325_tp+mx*( &
                  70621.44088156540229_tp+mx*( &
                  358977.2665825309926_tp+mx*( &
                  1.847238263723971684e6_tp+mx*( &
                  9.600515416049214109e6_tp+mx*( &
                  5.030767708502366879e7_tp+mx*( &
                  2.654441886527127967e8_tp+mx*( &
                  1.408862325028702687e9_tp+mx*( &
                  7.515687935373774627e9_tp)))))))))))))))
    else
            mx =  m-0.875_tp
            tp_elk=  2.473596173751343912_tp+mx*( &
                  3.727624244118099310_tp+mx*( &
                  15.60739303554930496_tp+mx*( &
                  84.12850842805887747_tp+mx*( &
                  506.9818197040613935_tp+mx*( &
                  3252.277058145123644_tp+mx*( &
                  21713.24241957434256_tp+mx*( &
                  149037.0451890932766_tp+mx*( &
                  1.043999331089990839e6_tp+mx*( &
                  7.427974817042038995e6_tp+mx*( &
                  5.350383967558661151e7_tp+mx*( &
                  3.892498869948708474e8_tp+mx*( &
                  2.855288351100810619e9_tp+mx*( &
                  2.109007703876684053e10_tp+mx*( &
                  1.566998339477902014e11_tp+mx*( &
                  1.170222242422439893e12_tp+mx*( &
                  8.777948323668937971e12_tp+mx*( &
                  6.610124275248495041e13_tp+mx*( &
                  4.994880537133887989e14_tp+mx*( &
                  3.785974339724029920e15_tp)))))))))))))))))))
      endif arg    
    return
  end function tp_elk
  
  
  pure function qp_elk( mc )
    real(qp),intent(in)      :: mc
    real(qp)                 :: qp_elk
    
    qp_elk = 0._qp
    return
  end function qp_elk  
  
  
  pure function dp_elk( mc )
    !Functions variables
    real(dp),intent(in)      :: mc
    real(dp)                 :: dp_elk
    !Working variables
    real(dp)                 :: m , mx , P , Q , kkc , nome 
    !Parameters
    real(dp),parameter       :: very_tiny = 1.e-99_dp
    real(dp),parameter :: D01 =                      real(1_i8,dp)/real(16_i8,dp)
    real(dp),parameter :: D02 =                      real(1_i8,dp)/real(32_i8,dp) 
    real(dp),parameter :: D03 =                     real(21_i8,dp)/real(1024_i8,dp)
    real(dp),parameter :: D04 =                     real(31_i8,dp)/real(2048_i8,dp) 
    real(dp),parameter :: D05 =                   real(6257_i8,dp)/real(524288_i8,dp)
    real(dp),parameter :: D06 =                  real(10293_i8,dp)/real(1048576_i8,dp) 
    real(dp),parameter :: D07 =                 real(279025_i8,dp)/real(33554432_i8,dp)
    real(dp),parameter :: D08 =                 real(483127_i8,dp)/real(67108864_i8,dp)
    real(dp),parameter :: D09 =              real(435506703_i8,dp)/real(68719476736_i8,dp)
    real(dp),parameter :: D10 =              real(776957575_i8,dp)/real(137438953472_i8,dp)
    real(dp),parameter :: D11 =            real(22417045555_i8,dp)/real(4398046511104_i8,dp)
    real(dp),parameter :: D12 =            real(40784671953_i8,dp)/real(8796093022208_i8,dp)
    real(dp),parameter :: D13 =          real(9569130097211_i8,dp)/real(2251799813685248_i8,dp)
    real(dp),parameter :: D14 =         real(17652604545791_i8,dp)/real(4503599627370496_i8,dp)
    !The following coefficients are not needed for double precision
    !real(dp),parameter :: D15 =        real(523910972020563_i8,dp)/real(144115188075855872_i8,dp)
    !real(dp),parameter :: D16 =        real(976501268709949_i8,dp)/real(288230376151711744_i8,dp)
    !!!Integer overflow for integer_kind=8 for the following coefficients
    !real(dp),parameter :: D17 =             935823746406530603._dp/295147905179352825856._dp
    !real(dp),parameter :: D18 =            1758220447807291611._dp/590295810358705651712._dp
    !real(dp),parameter :: D19 =           53030538453624441751._dp/18889465931478580854784._dp
    !real(dp),parameter :: D20 =          100268465197007602645._dp/37778931862957161709568._dp
         
    m = 1.0_dp - mc
arg:    if ( abs(m) < 1.e-16_dp ) then
          dp_elk = PIHALF
        elseif ( mc < very_tiny ) then
          dp_elk = LN4 - 0.5_dp*log(very_tiny)
        elseif ( mc < 1.11d-16 ) then
          dp_elk = LN4 - 0.5_dp*log(mc)
        elseif ( mc < 0.1_dp ) then
          nome = mc*(D01 + mc*(D02 + mc*(D03 + mc*(D04 + mc*(D05 + mc*(D06 &
               + mc*(D07 + mc*(D08 + mc*(D09 + mc*(D10 + mc*(D11 + mc*(D12 &
               + mc*(D13 + mc*D14)))))))))))))
          mx   = mc - 0.05_dp
!
!      K'
!
        kkc=1.591003453790792180_dp+mx*( &
            0.416000743991786912_dp+mx*( &
            0.245791514264103415_dp+mx*( &
            0.179481482914906162_dp+mx*( &
            0.144556057087555150_dp+mx*( &
            0.123200993312427711_dp+mx*( &
            0.108938811574293531_dp+mx*( &
            0.098853409871592910_dp+mx*( &
            0.091439629201749751_dp+mx*( &
            0.085842591595413900_dp+mx*( &
            0.081541118718303215_dp))))))))))
        dp_elk=-kkc*PIINV*log( nome )
    elseif(m <= 0.1_dp) then
            mx =  m-0.05_dp
            dp_elk = 1.591003453790792180_dp+mx*( &
                  0.416000743991786912_dp+mx*( &
                  0.245791514264103415_dp+mx*( &
                  0.179481482914906162_dp+mx*( &
                  0.144556057087555150_dp+mx*( &
                  0.123200993312427711_dp+mx*( &
                  0.108938811574293531_dp+mx*( &
                  0.098853409871592910_dp+mx*( &
                  0.091439629201749751_dp+mx*( &
                  0.085842591595413900_dp+mx*( &
                  0.081541118718303215_dp))))))))))
    elseif(m <= 0.2_dp) then
            mx  = m-0.15_dp
            dp_elk = 1.635256732264579992_dp+mx*( &
                  0.471190626148732291_dp+mx*( &
                  0.309728410831499587_dp+mx*( &
                  0.252208311773135699_dp+mx*( &
                  0.226725623219684650_dp+mx*( &
                  0.215774446729585976_dp+mx*( &
                  0.213108771877348910_dp+mx*( &
                  0.216029124605188282_dp+mx*( &
                  0.223255831633057896_dp+mx*( &
                  0.234180501294209925_dp+mx*( &
                  0.248557682972264071_dp+mx*( &
                  0.266363809892617521_dp+mx*( &
                  0.287728452156114668_dp))))))))))))
    elseif(m <= 0.3_dp) then
            mx =  m-0.25_dp
            dp_elk = 1.685750354812596043_dp+mx*( &
                  0.541731848613280329_dp+mx*( &
                  0.401524438390690257_dp+mx*( &
                  0.369642473420889090_dp+mx*( &
                  0.376060715354583645_dp+mx*( &
                  0.405235887085125919_dp+mx*( &
                  0.453294381753999079_dp+mx*( &
                  0.520518947651184205_dp+mx*( &
                  0.609426039204995055_dp+mx*( &
                  0.724263522282908870_dp+mx*( &
                  0.871013847709812357_dp+mx*( &
                  1.057652872753547036_dp)))))))))))
    elseif(m <= 0.4_dp) then
            mx =  m-0.35_dp
            dp_elk = 1.744350597225613243_dp+mx*( &
                  0.634864275371935304_dp+mx*( &
                  0.539842564164445538_dp+mx*( &
                  0.571892705193787391_dp+mx*( &
                  0.670295136265406100_dp+mx*( &
                  0.832586590010977199_dp+mx*( &
                  1.073857448247933265_dp+mx*( &
                  1.422091460675497751_dp+mx*( &
                  1.920387183402304829_dp+mx*( &
                  2.632552548331654201_dp+mx*( &
                  3.652109747319039160_dp+mx*( &
                  5.115867135558865806_dp+mx*( &
                  7.224080007363877411_dp))))))))))))
    elseif(m <= 0.5_dp) then
            mx  = m-0.45_dp
            dp_elk = 1.813883936816982644_dp+mx*( &
                  0.763163245700557246_dp+mx*( &
                  0.761928605321595831_dp+mx*( &
                  0.951074653668427927_dp+mx*( &
                  1.315180671703161215_dp+mx*( &
                  1.928560693477410941_dp+mx*( &
                  2.937509342531378755_dp+mx*( &
                  4.594894405442878062_dp+mx*( &
                  7.330071221881720772_dp+mx*( &
                  11.87151259742530180_dp+mx*( &
                  19.45851374822937738_dp+mx*( &
                  32.20638657246426863_dp+mx*( &
                  53.73749198700554656_dp+mx*( &
                  90.27388602940998849_dp)))))))))))))
    elseif(m <= 0.6_dp) then
            mx  = m-0.55_dp
            dp_elk = 1.898924910271553526_dp+mx*( &
                  0.950521794618244435_dp+mx*( &
                  1.151077589959015808_dp+mx*( &
                  1.750239106986300540_dp+mx*( &
                  2.952676812636875180_dp+mx*( &
                  5.285800396121450889_dp+mx*( &
                  9.832485716659979747_dp+mx*( &
                  18.78714868327559562_dp+mx*( &
                  36.61468615273698145_dp+mx*( &
                  72.45292395127771801_dp+mx*( &
                  145.1079577347069102_dp+mx*( &
                  293.4786396308497026_dp+mx*( &
                  598.3851815055010179_dp+mx*( &
                  1228.420013075863451_dp+mx*( &
                  2536.529755382764488_dp))))))))))))))
    elseif(m <= 0.7_dp) then
            mx =  m-0.65_dp
            dp_elk = 2.007598398424376302_dp+mx*( &
                  1.248457231212347337_dp+mx*( &
                  1.926234657076479729_dp+mx*( &
                  3.751289640087587680_dp+mx*( &
                  8.119944554932045802_dp+mx*( &
                  18.66572130873555361_dp+mx*( &
                  44.60392484291437063_dp+mx*( &
                  109.5092054309498377_dp+mx*( &
                  274.2779548232413480_dp+mx*( &
                  697.5598008606326163_dp+mx*( &
                  1795.716014500247129_dp+mx*( &
                  4668.381716790389910_dp+mx*( &
                  12235.76246813664335_dp+mx*( &
                  32290.17809718320818_dp+mx*( &
                  85713.07608195964685_dp+mx*( &
                  228672.1890493117096_dp+mx*( &
                  612757.2711915852774_dp))))))))))))))))
    elseif(m <= 0.8_dp) then
            mx =  m-0.75_dp
            dp_elk = 2.156515647499643235_dp+mx*( &
                  1.791805641849463243_dp+mx*( &
                  3.826751287465713147_dp+mx*( &
                  10.38672468363797208_dp+mx*( &
                  31.40331405468070290_dp+mx*( &
                  100.9237039498695416_dp+mx*( &
                  337.3268282632272897_dp+mx*( &
                  1158.707930567827917_dp+mx*( &
                  4060.990742193632092_dp+mx*( &
                  14454.00184034344795_dp+mx*( &
                  52076.66107599404803_dp+mx*( &
                  189493.6591462156887_dp+mx*( &
                  695184.5762413896145_dp+mx*( &
                  2.567994048255284686d6+mx*( &
                  9.541921966748386322d6+mx*( &
                  3.563492744218076174d7+mx*( &
                  1.336692984612040871d8+mx*( &
                  5.033521866866284541d8+mx*( &
                  1.901975729538660119d9+mx*( &
                  7.208915015330103756d9)))))))))))))))))))
    elseif(m <= 0.85_dp) then
            mx =  m-0.825_dp
            dp_elk=  2.318122621712510589_dp+mx*( &
                  2.616920150291232841_dp+mx*( &
                  7.897935075731355823_dp+mx*( &
                  30.50239715446672327_dp+mx*( &
                  131.4869365523528456_dp+mx*( &
                  602.9847637356491617_dp+mx*( &
                  2877.024617809972641_dp+mx*( &
                  14110.51991915180325_dp+mx*( &
                  70621.44088156540229_dp+mx*( &
                  358977.2665825309926_dp+mx*( &
                  1.847238263723971684e6_dp+mx*( &
                  9.600515416049214109e6_dp+mx*( &
                  5.030767708502366879e7_dp+mx*( &
                  2.654441886527127967e8_dp+mx*( &
                  1.408862325028702687e9_dp+mx*( &
                  7.515687935373774627e9_dp)))))))))))))))
    else
            mx =  m-0.875_dp
            dp_elk=  2.473596173751343912_dp+mx*( &
                  3.727624244118099310_dp+mx*( &
                  15.60739303554930496_dp+mx*( &
                  84.12850842805887747_dp+mx*( &
                  506.9818197040613935_dp+mx*( &
                  3252.277058145123644_dp+mx*( &
                  21713.24241957434256_dp+mx*( &
                  149037.0451890932766_dp+mx*( &
                  1.043999331089990839e6_dp+mx*( &
                  7.427974817042038995e6_dp+mx*( &
                  5.350383967558661151e7_dp+mx*( &
                  3.892498869948708474e8_dp+mx*( &
                  2.855288351100810619e9_dp+mx*( &
                  2.109007703876684053e10_dp+mx*( &
                  1.566998339477902014e11_dp+mx*( &
                  1.170222242422439893e12_dp+mx*( &
                  8.777948323668937971e12_dp+mx*( &
                  6.610124275248495041e13_dp+mx*( &
                  4.994880537133887989e14_dp+mx*( &
                  3.785974339724029920e15_dp)))))))))))))))))))
      endif arg
    return 
  end function dp_elk
end module complete_elliptic
